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1  Abstract 


This  paper  investigates  the  validity  of  commonly  used  terramechanics  models  for  light-weight  vehicle  applications 
while  accounting  for  experimental  variability.  This  is  accomplished  by  means  of  cascading  uncertainty  up  to  the 
terminal  point  of  operations  measurement.  Vehicle -terrain  interaction  is  extremely  complex,  and  thus  models  and 
simulation  methods  for  vehicle  mobility  prediction  are  largely  based  on  empirical  test  data.  Analytical  methods  are 
compared  to  experimental  measurements  of  key  operational  parameters  such  as  drawbar  force,  torque,  and  sinkage. 
Models  of  these  operational  parameters  ultimately  depend  on  a  small  set  of  empirically  determined  soil  parameters, 
each  with  an  inherent  uncertainty  due  to  test  variability.  The  soil  parameters  associated  with  normal  loads  are 
determined  by  fitting  the  dimensionless  form  of  Bekker’s  equation  to  the  data  given  by  the  pres  sure -sinkage  test.  In 
a  similar  approach,  the  soil  parameters  associated  with  shear  loads  are  determined  by  fitting  Janosi  and  Hanamoto’s 
equation  to  the  data  given  by  the  direct  shear  test.  An  uncertainty  model  is  used  to  propagate  the  soil  parameter 
variability  through  to  the  wheel  performance  based  on  Wong  and  Reece.  The  commonly  used  analytical  model  is 
shown  to  be  inaccurate  as  the  envelope  of  model  uncertainty  does  not  lie  within  the  experimental  measures, 
suggesting  that  model  improvements  are  required  to  accurately  predict  the  performance  of  light-weight  vehicles  on 
deformable  terrain. 


2  Keywords 


Terramechanics;  Robotic  vehicles;  Bekker  and  Wong  models;  Soil  test  bed;  Uncertainty  propagation;  Stochastic 
modeling 


3  Background 


The  study  of  the  interaction  of  wheeled  and  tracked  vehicles  with  natural  terrain  is  dominated  by  the  discipline  of 
terramechanics.  Terramechanics  research  over  the  past  50  years  has  primarily  focused  on  large,  heavy  military 
vehicles.  A  substantial  body  of  terramechanics  research  has  been  performed  at  the  U.S.  Army  Tank  Automotive 
Research,  Development,  and  Engineering  Center  (TARDEC)  and  the  U.S.  Army  Engineer  Research  and 
Development  Center  (ERDC)  that  led  to  the  development  of  various  mobility  prediction  methodologies  including 
the  NATO  Reference  Mobility  Model  (NRMM).  These  methodologies  are  numerical  algorithms  for  predicting 
cross-country  vehicle  movement  at  length  scales  of  several  meters  to  several  kilometers.  They  are  based  on 
empirical  results  drawn  from  years  of  resource-intensive  experimental  testing  and  have  been  used  widely  by  the 
military  community.  However,  as  a  consequence  of  their  empirical  nature,  while  the  methods  are  useful  for 
prediction  of  large,  heavy  vehicle  mobility,  it  remains  an  open  question  whether  they  can  be  reliably  used  for  the 
prediction  of  small,  lightweight  vehicle  mobility. 

Recently  the  Department  of  Defense  (DoD)  has  devoted  substantial  resources  toward  the  development  of  small, 
lightweight  ground  vehicles.  These  vehicles  are  often  less  than  36  inches  in  length  and  weigh  less  than  100  lb.  They 
are  equipped  with  wheels,  tracks,  or  bio-inspired  limb-like  appendages.  Due  to  the  lack  of  mobility  models  for  these 
vehicles,  there  is  a  lack  of  simulation  methods  for  this  class  of  vehicles.  Also,  since  these  vehicles  have  only 
recently  been  adopted  by  the  DoD,  there  is  a  lack  of  systematic  empirical  test  data.  As  a  result,  these  vehicles  have 
primarily  been  developed  based  on  ad  hoc  design  rules,  limited  empirical  testing,  and  application  of  classical  Bekker 
theory.  Thus,  there  is  a  lack  of  methods  (based  on  simulation  or  analysis)  to  reliably  predict  the  mobility  and 
performance  of  these  systems. 

This  paper  will  describe  two  experimental  methods  used  to  characterize  mechanical  soil  behavior  for  lightweight 
vehicles.  A  pressure-sinkage  test  and  a  standard  direct  shear  test  as  outlined  by  Bekker  in  [1]  were  performed  on  a 
cohesion-less  soil  [2].  This  paper  investigates  identification  of  soil  parameters  from  the  experimental  data  and 
assesses  the  data  variability,  which  is  inherent  in  both  tests.  The  variability  in  the  soil  parameters  is  propagated  to 
determine  the  overall  wheel  performance  uncertainty.  Probabilistic,  rather  than  strictly  deterministic,  soil  behavior  is 
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considered  in  this  paper  because  it  is  an  increasingly  important  consideration  for  accurate  modeling  of  lightweight 
robotic  vehicle  performance. 


4  Methodology 


Semi-empirical  methods  for  modeling  wheel  performance,  like  the  one  used  in  this  paper  by  Wong  and  Reece  [3], 
rely  on  the  relationship  between  soil  sinkage  and  resistance  force  to  infer  the  normal  stress  under  a  wheel.  To  predict 
the  tractive  force,  the  shearing  strength  of  the  soil  is  analyzed  based  on  Coulomb’s  formula  [4].  Methods  of  this 
class  are  ultimately  based  on  experimentally  determined  soil  parameters,  whose  inherent  variability  causes 
uncertainty  in  the  determination  of  wheel  performance.  Section  4.1  describes  the  equipment  and  methods  used  to 
determine  the  variability  in  soil  parameters;  Section  4.2  describes  the  wheel  model  in  which  these  parameters  are 
used  and  the  techniques  used  to  compare  prediction  to  experimental  data. 

4. 1  Soil  Parameter  Variability 

Mojave  Martian  Simulant  (MMS)  [2]  was  employed  as  a  test  medium  for  the  experiments  in  this  paper.  MMS  is  a 
mixture  of  finely  crushed  and  sorted  granular  basalt  intended  to  mimic,  both  at  a  chemical  and  mechanical  level, 
Mars  soil  characteristics.  No  direct  application  to  Mars  rovers  is  provided  in  this  paper,  the  simulant  is  a  frictional 
soil  found  in  dry,  sandy  terrain.  The  particle  size  distribution  of  MMS  spans  from  the  micron  to  millimeter  level 
with  80%  of  particles  above  the  10  micron  threshold. 

4.1.1  Experimental  Equipment 

4.1.1.1  Pressure-Sinkage  Test 

The  sinkage  characteristics  of  the  MMS  were  measured  using  the  pressure-sinkage  test  shown  in  Figure  1.  The 
pressure-sinkage  test  used  a  plate  to  penetrate  the  soil  under  controlled  test  conditions,  while  pressure  and  sinkage  of 
the  plate  were  directly  measured.  A  series  of  tests  and  various  plate  sizes  allowed  an  investigation  of  both  the 
influence  of  the  pressure-sinkage  parameters  as  described  in  Wong’s  methodology  [5]  and  test-to-test  variability. 

The  test  unit  was  designed  to  systematically  penetrate  the  soil  with  a  downward  velocity  of  10  mm/s.  Penetration 
tests  were  performed  with  three  different-sized  rectangular  plates  (3x15,  5x15,  and  7x15  cm2).  The  tests  were 
repeated  15  times  for  each  plate.  A  load  cell  and  draw-wire  encoder  recorded  the  force  and  corresponding  sinkage 
during  each  test.  Between  tests,  the  soil  was  loosened  with  a  stick  and  then  leveled  to  return  the  MMS  to  a  nominal 
density  of  1.7  g/cm3.  Figure  2  shows  a  picture  of  a  penetration  plate  about  to  enter  the  soil. 


Figure  1:  CAD  drawing  of  the  pressure-sinkage  test  rig.  A  Figure  2:  Penetration  plate  and  soil  (side  view). 

load  cell  and  draw  wire  encoder  recorded  the  force  and 
sinkage  of  a  plate  that  was  pressed  into  the  soil  at  10  mm/s. 
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4.1.1.2  Direct  Shear  Test 

The  direct  shear  test,  shown  in  Figure  3,  is  used  to  measure  the  shear  strength  properties  of  the  MMS,  specifically 
the  cohesion,  angle  of  friction,  and  shear  modulus.  A  sample  of  the  soil  is  contained  in  between  two  rigid  discs  that 
are  held  in  place  by  a  shear  box.  The  shear  box  is  aligned  under  a  load  cell  that  applies  a  normal  force  to  the  soil. 


loading  cap 


porous 

plate 


upper  box 


lower  box 
soil  sample 

porous  plate 


Figure  3:  Experimental  device  for  performing  the  direct  shear  test. 

The  load  cell  is  attached  to  a  vertical  translational  joint  that  uses  a  linear  variable  differential  transformer  to  measure 
displacement  of  the  soil.  The  top  of  the  shear  box  is  clamped  so  that  the  lower  half  can  be  moved.  The  horizontal 
force  required  to  displace  the  soil  horizontally  is  measured  by  a  dynamometer.  The  applied  vertical  force  and 
measured  horizontal  force  can  be  transformed  into  the  normal  and  shear  stress,  respectively.  The  horizontal  and 
vertical  soil  displacement  is  also  output. 

4.1.2  Parameter  Estimation 


4.1.2.1  Pressure-Sinkage  Parameter  Estimation 

To  describe  the  pressure-sinkage  relationship  of  a  soil’s  deformation  under  a  rectangular  plate,  Bekker  suggested 
Equation  (1)  based  on  the  soil  mechanics  work  originally  performed  by  Terzaghi  [6]: 


K  7 
+K 

b  9 


z  , 


(1) 


where  a  is  pressure,  b  is  plate  width,  z  is  soil  sinkage,  and  n  ,  kc  and  k  ,  are  soil  parameters.  The  three 

parameters  are  empirical  and  have  no  intrinsic  physical  meaning.  The  parameters  may  be  estimated  from 
experimental  pressure-sinkage  data  if  at  least  two  different  plate  sizes  are  tested. 

The  process  of  determining  the  soil  parameters  from  experimental  data  was  originally  performed  by  Wong  using  a 
weighted  least  squares  method  [7],  which  is  the  primary  estimation  technique  used  in  this  paper.  From  Equation  (1), 
taking  the  logarithms  of  both  sides  gives: 


In  <j  -  In 


—  +  ka 
b  9 


+  n  In  Z- 


(2) 


An  error  function  is  defined  as: 
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U  'J  J 

where  wr  is  a  weighting  factor.  Standard  least  squares  minimization  of  equation  (3)  would  set  wr  equal  to  1. 
However,  as  Wong  explains,  error  would  be  biased  toward  low  pressures  since  the  actual  measured  values  are  cr 
and  z  ,  not  ln  cr  and  lnz  .  To  account  for  minimization  in  the  log-log  domain  and  give  equal  weight  to  all  data 
points,  Wong  sets  the  weighting  factor  wr=cr2 . 

Since  kc,  b  ,  and  are  constant  values  for  a  given  test,  they  can  be  replaced  by  a  single  constant: 


(4) 


Wong  provides  the  following  cost  function  over  N  data  points  that  is  minimized  to  find  optimal  n  and  keq  in 
Equation  (5): 

F  =  Z ai  [ln ai  “ln(/c«,)-»lns]2-  (5) 

i=\ 

The  optimality  condition  leads  to  the  following  two  equations: 

=  -2^^cr2  ln  cr  ln  z  -  In  keq  y<T2  ln  z-n^jr2  (lnz)2  J  =  0  (6) 


dF 

Kq 


—  [JV  In  cr + 2  ln  keq  +  2 ri^jj2  ln  z ]  =  0. 


(7) 


Solving  equations  (6)  and  (7)  for  n  and  keq ,  yields: 

y/T2  »2  ln  cr  In  z  -  ^cr2  ln  cr  y<r2  ln  z 

Z-2Z-2(m^-&2lnz)2 


y<r2  In (T-n'Yp'2  ln ; 


If  the  above  analysis  is  performed  for  two  sets  of  data  corresponding  to  two  different  plate  sizes 
(blf  b2 ),  the  average  n  and  two  corresponding  values  of  keq  can  be  found: 

(bL+("L 


(8) 


(9) 


(10) 
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(b  )  - 
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(12) 


Using  the  two  values  of  keq ,  kc  and  kip  can  be  derived: 


U2  u\ 


(13) 


k  =  (^L,  -(UJ- 


(14) 


Wong’s  methodology  is  a  clear  and  succinct  way  to  obtain  the  three  pressure-sinkage  soil  parameters  (n  ,kc,k  ) 
from  a  set  of  experimental  data. 

To  avoid  difficulties  resulting  from  covariance  between  soil  parameters  (shown  by  example  in  Section  5. 1.1.1), 
Reece’s  revision  of  Bekker’s  equation  [3]  was  implemented  by  the  authors  of  this  paper.  Reece’s  equation  is: 


a  = 


(15) 


where  kx  and  k2  are  empirical  soil  parameters.  In  simplified  terms,  Equation  (15)  reduces  to: 


(16) 


where  the  k  coefficient  is  denoted  as  keq  to  distinguish  from  the  keq  of  Bekker’s  equation.  To  provide  a  more 

thorough  analysis  of  pressure-sinkage  parameter  estimation  methods,  four  additional,  modified  approaches  were  also 
considered.  The  resulting  parameters  of  all  five  methods  are  compared.  Because  Wong’s  pressure-sinkage  relation 
and  Reece’s  modification  have  the  same  format  for  any  single  plate,  the  results  of  the  following  methods  apply 
equally  to  both. 

Method  1  is  performed  using  Wong’s  equations  (8)-(14)  previously  given  for  individual  experimental  tests. 

Method  2  is  a  modified  version  of  the  technique  presented  by  Wong.  The  weighting  factor  wr  in  equation  (3)  is 
changed  from  cr2  to  1 .  The  resulting  equations  for  n  and  keq  are  derived  analytically  using  a  similar  approach  as 
Wong: 


n  = 


^  In  cr  In  z  - 

S(lnz)2 


N 


(IH 

N 


(17) 
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Kq  =eXP 


^Incr-n^ln: 

N 


(18) 


From  equation  (17),  kc  and  can  be  derived  in  the  same  manner  as  equations  (13)  and  (14). 

For  the  final  three  methods,  Levenberg-Marquardt  numerical  optimization  (from  the  lsqnonlin  ()  function  in 
MATLAB’s  optimization  suite  [8])  is  used  to  find  n ,  keq  pairs  that  minimize  the  respective  error  functions. 

Method  3  minimizes  error  in  the  ordinary  z  —  cr  domain.  No  weighting  was  considered  because  error  and  actual 
measurements  are  in  the  same  domains.  The  minimization  function  is: 

F=ZU-V"]2-  (19> 


Method  4  has  an  identical  error  function  to  that  of  method  2:  both  are  in  the  log-log  domain  and  have  wr  set  to  1. 
The  only  difference  is  that  method  4  is  numerically  minimized,  whereas  method  2  is  analytically  derived.  The 
minimization  function  is: 


F=J^[huj-lnkx/ 


(20) 


Method  5  is  the  numerical  version  of  Wong’s  method,  method  1.  The  minimization  function  is: 

F  =^jr2  [lncr-ln/c^  -n\nz~\  . 


(21) 


Methods  1-5  are  summarized  in  Table  1. 

Table  1:  Alternative  methods  to  estimate  pressure-sinkage  parameters. 


Method 

Type 

Domain 

Weighting  factor,  wr 

1.  Wong’s 

Analytical 

log-log 

<72 

2.  Wong’s  (modified) 

Analytical 

log-log 

1 

3.  Least-square-curve-fit 

Numerical 

Z-CT 

N/A  (z  —  cr domain) 

4.  Least-square-curve-fit 

Numerical 

log-log 

1 

5.  Least-square-curve-fit 

Numerical 

log-log 

cr2 

Individual  parameter  estimation  of  the  45  pressure-sinkage  tests  was  performed  using  each  of  the  five  methods.  The 
error  of  a  fit  for  all  methods  is  calculated  using  Wong’s  suggestion: 


J  N- 2 


N 


(22) 
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where  <7m  is  measured  pressure,  oc  is  estimated  pressure,  and  N  is  the  number  of  data  points. 

4.1.2. 2  Direct  Shear  Parameter  Estimation 

Several  soil  parameters  can  be  determined  using  the  direct  shear  test.  Specifically,  the  residual  shear  stress,  zres ,  and 
shear  modulus,  K  ,  can  be  determined  by  fitting  the  Janosi  and  Hanamoto  equation  to  each  data  set,  and  the 
cohesion,  c  ,  angle  of  internal  friction,  <p  ,  can  be  found  using  the  Mohr-Coulomb  failure  criteria. 

4.1. 2. 2.1  Janosi-Hanamoto  Equation  for  Shear  Stress 

The  shear-displacement  expression  suggested  by  Janosi  and  Hanamoto  [9]  assumed  the  form  of  Equation  (23)  for 
loose  soils: 


t  =  t„ 


1-e1 


(23) 


where  z  represents  the  shear  stress  and  j  is  the  displacement  due  to  shearing.  The  parameters  zres  and  K  are 
determined  by  minimizing  the  sum  of  the  differences  between  the  experimental  value  of  r  and  the  estimate  at  the 
ith  data  point  in  Equation  (24)  using  the  Levenberg-Marquardt  algorithm: 


)=I 


f 

f  ~ii  > 

\ 

z.-z 

l-eK 

i  res 

V 

V  J 

/ 

(24) 


The  optimality  condition  leads  to  the  following  two  equations: 


(25) 


(26) 


where  the  partial  derivatives  are  determined  using  a  numerical  central  differences  approach. 

4.1. 2. 2. 2  Mohr-Coulomb  Failure  Criteria 

For  a  given  normal  load,  the  soil  is  said  to  fail  when  it  reaches  its  residual  shear  stress.  According  to  the  Mohr- 
Coulomb  failure  criteria  [10],  a  line  of  best  fit  can  be  determined  by  plotting  the  residual  shear  stress,  zres ,  as  a 
function  of  normal  stress  ,  a  ,  as  follows: 


Tres  =c  +  cr  tan  <p. 


(27) 


where  the  intercept  is  equal  to  the  cohesion  and  the  slope  is  related  to  the  angle  of  internal  friction.  The  line  of  best 
fit  is  determined  using  the  closed-form  equation  for  linear  fitting  over  N  data  points  (Equation  (28)).  The  closed- 
form  equation  will  be  useful  for  propagating  the  uncertainty  in  Section  4.2.3: 
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N 

i= 1 

N 

If. 

i= i 

- 1 

cj  p 

-s 

N 

V  G  T 

i  res,i 

i= 1 

N 

N 

fcf 

w? 

_ 1 

N 

Vr 

res,i 

i= 1 

(28) 


4.2  Wheel  Performance  Uncertainty  Propagation  and  Verification 

This  section  details  the  experimental  method  used  to  obtain  the  wheel  performance  measurements,  explains  the 
theory  behind  the  prediction  of  wheel  performance,  and  describes  the  derivation  behind  the  method  of  uncertainty 
propagation.  Since  the  techniques  used  in  this  paper  rely  on  estimates  of  the  stress  distribution  to  determine  wheel 
performance,  the  ability  to  measure  these  quantities  allow  for  a  direct  comparison  of  prediction  to  reality.  The 
uncertainty  propagation  gives  an  envelope  of  system  performance  that  the  experimental  data  can  be  compared 
against. 

4.2.1  Experimental  Methodology 

A  complete  investigation  into  the  performance  of  the  wheel  requires  insight  into  several  different  measures  of  the 
wheel-soil  interaction.  In  this  paper,  physical  measurements  of  contact  stresses,  force,  torque,  and  sinkage  were  used 
to  compare  to  the  theoretical  model.  A  single-wheel  test  rig,  shown  in  Figure  4,  was  used  to  empirically  investigate 
the  wheel  motion  under  controlled  wheel  slip  and  normal  loading  conditions  on  the  cohesion-less  soil  [11].  This  test 
rig  enables  the  control  of  velocities  and  application  of  loads  through  interchangeable  running  gear  within  a  confined 
soil  bin  of  dimensions  1.5  m  long,  0.7  m  wide,  and  0.4  m  deep.  The  drawbar  pull,  wheel  torque,  and  sinkage  were 
measured  for  a  lug-less  rigid  wheel  with  a  radius  of  0.13  m  and  a  width  of  0.16  m  at  slip  ratios  varying  from  zero  to 
unity.  Tests  were  run  for  two  different  cases:  normal  loads  of  70  N  and  135  N. 


Figure  4:  A  single  wheel  test  bed  was  used  to  measure  the  drawbar  pull,  torque,  and  sinkage  at  varying  slips  and  normal  loads. 
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Figure  5:  A  wheel  with  a  custom  force  sensor  array  was  used  to  measure  the  stress  distribution  at  the  contact  patch. 


To  measure  the  stress  distribution  at  the  contact  patch  of  the  wheel,  an  identical  wheel  with  a  custom  force  sensor 
array  located  at  the  wheel  surface  was  used  [12].  The  force  sensors  are  strain  gauge -based  flexural  elements  with 
interchangeable  interface  surfaces  that  are  designed  for  integration  with  wheels  or  other  running  gear.  The  normal 
and  shear  stresses  are  estimated  based  on  5  sensors,  shown  in  Figure  6,  that  are  located  at  discrete  points  from  the 
center  to  the  edge  of  the  wheel-soil  interface. 


Wheel  top  view 


1  Sensor  1  Sensor 

1  Sensor  1 

1  Sensor  1 

1  1  1" 

[wj 

[jvj 

11 

Figure  6:  A  top  view  of  the  wheel  shows  the  layout  of  the  strain  gauge-based  sensors.  Note  that  the  sensors  are  arranged  with 
Sensor  I  at  the  center  of  the  wheel  and  Sensor  V  at  the  outer  edge. 

Each  sensor,  shown  in  Figure  7,  is  a  solid-state  L-shaped  aluminum  flexure  instrumented  with  two  full  bridge  strain 
gages.  The  sensor  is  mounted  rigidly  to  the  running  gear,  and  its  interface  element  is  exposed  to  the  soil.  The 
interface  element  is  generally  subjected  to  normal,  N  ,  and  shear,  T  ,  loading.  These  forces  cause  the  flexure 
elements  to  deflect  in  a  linear  elastic  manner.  From  measured  deflection,  and  given  prior  calibration  data,  the 
applied  forces  can  be  uniquely  computed,  with  axial  strain  intrinsically  rejected  by  the  full  bridge  configuration. 
Stress  can  then  be  inferred  assuming  uniform  pressure  distribution  over  the  known  sensors’  head  area. 


T  cLx  +  N  d2 


Td± 

Figure  7:  Working  scheme  of  the  custom  force  sensor  for  interfacial  stress  measurement. 

4.2.2  Mathematical  Model  of  Rigid  Wheel  on  Soft  Soil 

The  rigid  wheel  free-body  diagram  based  on  the  work  of  Wong  and  Reece  [3],  shown  in  Figure  8,  is  used  in  this 
paper  to  model  the  interaction  between  the  wheel  and  the  soil.  Using  this  model,  the  drawbar  pull  D  ,  torque  T  ,  and 
sinkage  z  ,  can  be  estimated  for  a  wheel  of  weight  W  ,  radius  r ,  and  wheel  width  b  ,  travelling  at  a  linear  velocity 
v . 
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Figure  8:  Forces ,  torque ,  and  stresses  on  a  driven  rigid  wheel. 

The  sinkage  of  the  wheel  is  commonly  converted  into  polar  coordinates  using  the  wheel  hub  as  the  origin.  Once  the 
limits  of  the  contact  patch,  6X  and  02 ,  between  the  wheel  and  the  soil  are  determined,  the  drawbar  pull  and  torque 
can  be  calculated  by  integrating  the  radial  and  tangential  stresses  over  the  wheel. 

At  a  specified  slip,  i ,  it  is  possible  to  determine  the  sinkage,  forces,  and  torques  that  act  on  a  wheel  given  the  soil 
parameters  and  wheel  properties  [3].  A  force  balance  in  the  vertical  direction  yields  an  equation  for  the  weight,  W  , 
of  the  tire: 


W  =  rb  J  (crcos#  +  rsin#)<f$, 

-02 


(29) 


where  the  normal  stress,  a  ,  is  determined  by  Equation  (16)  and  the  shear  stress,  r  ,  is  computed  as: 

t  =  (c  +  cj  tan  cp )  (l  -  e~]IK  ) 


where  the  shear  deformation,  j,  is  based  on  the  velocity  of  the  slip: 

/  =  r[(#!  -#)-(!-/) (sin^  -sin  (9)] 


(30) 


(31) 


Based  on  geometry,  the  sinkage,  z  ,  can  be  converted  to  polar  coordinates  using: 

z  =  (cos0-cos0l)r  0M  <9<9l 


(32) 


However,  it  is  important  to  note  that  Equation  (32)  can  only  be  used  when  the  normal  stress  on  the  wheel  is 
increasing,  which  occurs  from  the  front  contact  angle,  9X  ,  to  a  maximum  radial  stress,  0M  .  At  0M  the  stress  will 

begin  to  decrease  in  a  similar  fashion  until  reaching  zero  at  the  rear  contact  angle,  02 .  The  symmetry  between  the 
front  and  rear  regions  is  employed  to  derive  a  relationship  for  sinkage  when  the  normal  stress  is  decreasing: 


z  = 


cos 

V  v 


01- 


6>  +  <92  ^ 

\0M  +  d2  j 


(0,-0*) 


-  cos  6X 


-e2<e<eu 


(33) 
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where  0M  can  be  approximated  as  a  function  of  0X  ,  i ,  and  two  coefficients  that  give  an  estimation  of  the  maximum 
stress,  cx  and  c2  [3]: 

oM={c i+'VM  (34> 

With  the  assumption  that  02  -  0  (which  is  true  when  rut  recovery  is  small),  the  only  unknown  in  Equation  (29) 
becomes  6X .  The  integral  in  Equation  (29)  cannot  be  solved  analytically  and  is  solved  numerically  via  an  iterative 

technique.  After  the  front  contact  angle  is  determined,  the  forces  (thrust  H  ,  motion  resistance  R  ,  and  drawbar  pull 
D  )  and  the  input  torque,  T  ,  can  be  determined  by  integrating  the  stresses  over  the  wheel  contact  area: 

H  =  rb  cos  OdO  (35) 

-02 

01 

R  =  rb  J  <j  sin  OdO  (36) 

-02 

D  =  H-R  (37) 

0i 

T  =  r2b  I  rd6  (38) 

-02 


Finally,  the  maximum  wheel  sinkage,  z0 ,  can  be  determined  using  Equation  (32)  at  0  =  0: 

Z0  =(l-cos  0x)r 


(39) 


4.2.3  Uncertainty  Propagation 

The  objective  of  any  uncertainty  analysis  is  to  obtain  an  estimate  of  the  overall  uncertainty,  i.e.  uncertainty 
envelope,  of  an  output  given  the  individual  uncertainties  of  the  inputs.  Previous  work,  including  Virginia  Tech’s 
Center  for  Vehicle  Systems  and  Safety  [13]  and  the  Massachusetts  Institute  of  Technology’s  Robotic  Mobility 
Group  [14],  used  Polynomial  Chaos  Theory  to  quantify  the  uncertainty  associated  with  terramechanics  systems.  The 
method  for  propagating  the  uncertainty  in  this  paper  is  described  in  [15],  and  summarized  as  follows: 

Given  an  arbitrary  function,  r  ,  that  depends  on  the  measured  variables  xx,x2,...  whose  distributions  are  normal: 

r  =  f(xl,x2,..),  (40) 


the  variance  of  r ,  cr  ,  can  be  determined  statistically  based  on  discrete  values  of  the  function,  ,  and  the  mean  of 
the  function,  r' ,  using: 


(41) 
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The  difference,  8rt ,  between  a  particular  value,  r  ,  and  its  mean  value,  r' ,  can  be  determined  from  a  Taylor  series 
expansion,  where  higher  order  terms  are  assumed  to  be  negligible: 


^  dr  ^ 


Sri  ={ri~r)= (*»  -*2) 


^  dr  ^ 


+  ... 


(42) 


If  the  equation  for  Sr  is  substituted  into  Equation  (41),  we  obtain: 


1  N 

<j)  =  lim  —  V 

n->  oo  yy 

1 Y  «=i 


(  '  \ 

fsk 

/  • \(  dr  ^ 

K  ~xi ) 

l 
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1 

£ 

+ 

+  ... 

(43) 


which  can  be  rearranged  to  yield: 


t  ■  \2 1 

(  dr)  ,  ,n2| 

f  dr  ) 

1 3*2  J 

r  N  AT 

1  y  1- 1 

/  -  \  /  .  \ 

r  dr  ^ 

r  dr  ^ 

+2\xu  -Xj  Jyx2i -x2) 

[dx2) 

(44) 


Note  that  the  variances  for  x,  and  x2  can  be  computed  as: 


cr  =  lim 

Xj  A^— >qo 


1  N  0 

-YU-x'f 

Jl  ]) 


(45) 


where  Equation  (45)  can  be  used  to  simplify  Equation  (44),  resulting  in  Equation  (46): 


__2  ^  2 

(7  =  (7 

r 


dr) 

2 

+  <Tx 

f  dr) 

2 

+  2cr  I 

f5rl 

f5rl 

U*iJ 

*2  | 

Vdx2) 

x\x2  | 

dxi) 

(46) 


Note  that  the  covariance  term,  <7x^  ,  is  typically  neglected  because  it  is  assumed  that  the  variables  x1  and  x2  are 

statistically  independent  of  one  another.  The  above  equation  gives  the  standard  deviation  of  the  independent 
variable  given  the  standard  deviation  of  the  variables  that  it  depends  on.  A  similar  statement  can  be  made  for  the 
uncertainty,  to  get  the  response  uncertainty: 


2  2 

11.  =  11 


dr  ) 

2 

2 

+  UZ 

r  dr  ^ 

2 

+  2  M 

r  dr  ^ 

r  dr  ^ 

kxj 

x2 

)dx2) 

xlx2  | 

[dxj 

{dx2j 

(47) 


Equation  (47),  hence,  defines  how  the  uncertainties  of  underlying  variables  propagates  into  the  overall  response.  The 
specific  parameters  that  contribute  to  the  uncertainty  of  the  wheel  performance  are  listed  in  Table  5  and  Table  6. 
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5  Results 


5.1  Soil  Parameter  Variability 

5.1.1  Pressure-Sinkage  Test 

A  total  of  45  pressure-sinkage  tests  were  performed  to  analyze  the  bearing  characteristics  of  the  soil.  The  tests  were 
repeated  15  times  for  three  rectangular  plates  with  different  widths,  b  (3,  5,  and  7  cm).  Parameters  were 
determined  using  the  five  different  methods  outlined  previously  in  Section  4. 1.2.1. 

Figure  9  shows  the  fits  produced  by  the  first  three  methods  for  a  single  test.  Methods  1  and  3  yield  similar  curves, 
though  varying  slightly.  The  curve  produced  by  method  2  appears  very  different  and  is  strongly  biased  towards 
fitting  the  initial  data  points  of  each  test.  Methods  1  and  2  produced  results  identical  to  methods  5  and  4, 
respectively.  The  same  trends  of  each  different  method  fit  are  seen  in  all  of  the  other  pressure-sinkage  tests. 


0  0.01  0.02  0.03  0.04 

Sinkage  [m] 

Figure  9:  Sample  estimation  for  5 -cm  plate. 

It  is  important  to  recognize  that  no  exponential  function  will  closely  fit  the  actual  shape  of  the  data  curves.  This  is  a 
direct  consequence  of  the  fact  that  Bekker’s  exponential  relationship  in  equation  (1)  is  an  empirical  estimation,  not  a 
physical  law,  and  it  does  not  account  well  for  deviance  in  the  early  behavior  (when  sinkage  is  less  than  about  15 
mm). 

Average  values  of  n  ,  k  ,  and  least-squares  error  for  each  fitting  method  and  plate  size  are  shown  in  Table  2. 


Table  2:  Average  parameters  and  error  for  3 -cm,  5-cm,  and  7 -cm  plates,  Bekker’s  equation. 


Plate  Width, 
b  (cm) 

3 

5 

7 

Method  # 

n 

^eq,3 

error 

n 

^eq,5 

error 

n 

keq,7 

error 

[MPa/mn] 

[MPa/mn] 

[MPa/mn] 

1,5 

1.01 

2.48 

0.040 

0.84 

1.90 

0.071 

1.082 

0.231 

0.096 

2,4 

1.39 

9.89 

0.196 

1.44 

25.7 

0.276 

1.114 

0.266 

0.220 

3 

1.03 

2.66 

0.039 

0.89 

2.21 

0.067 

0.976 

0.252 

0.092 

Page  14  of  33 


UNCLASSIFIED 


Methods  1,  3  and  5  have  relatively  small  fit  errors.  The  parameters  found  by  methods  2  and  4  are  substantially 
different  than  those  of  the  other  methods;  the  sinkage  exponent  n  is  much  higher  and  the  coefficient  keq is  much 
larger.  As  noted  previously,  these  methods  fit  low- sinkage  data  points  well  but  poorly  match  the  rest  of  the  data. 

Several  important  conclusions  may  be  drawn  from  these  results.  First,  Wong’s  method  (method  1)  has  proven  to  be 
very  reliable  as  compared  to  a  similar  numerical  method  (method  3).  Secondly,  it  gives  a  very  good  approximation 
to  the  best  fit  in  the  z  —  cr  domain.  Nevertheless,  estimation  in  the  log  domain  and  in  the  z  —  cr  domain  will  yield 
different  parameters.  A  third  observation  is  that  the  weighting  factor  a2  is  necessary  for  best  results  in  the  log-log 
domain. 

The  error  found  by  methods  1  and  3  is  comparable.  The  ( n ,  k  )  parameter  pairs  found  by  methods  1  and  3  also 

produce  a  similar  pressure-sinkage  curve,  so  it  is  concluded  that  either  method  may  be  used  for  comparable 
predictions.  It  is  interesting  to  note,  however,  that  the  individual  parameters  yielded  from  methods  1  and  3  can  be 
significantly  different  despite  both  yielding  a  good  fit. 

5. 1.1.1  Parameter  Covariance 

The  covariance  demonstrated  by  n  and  keq  requires  attention  before  proceeding  to  a  statistical  analysis.  Figure  10 

illustrates  the  apparent  correlation  of  Bekker’s  parameters  from  Section  4. 1.2.1.  Theoretically,  every  test  should 
yield  a  similar  n ,  but  each  plate  size  should  yield  a  different  keq. 


Figure  10:  Relationship  between  n  and  keq,  Bekker’s  equation. 

In  this  case,  n  and  keq  show  positive  correspondence.  This  correspondence  may  be  explained  by  noting  that 
modifying  the  units  of  keq  has  a  significant  effect.  The  relationship  between  n  and  keq  estimations  exist  because 
different  combinations  may  yield  a  similar  pressure  prediction,  and  because  the  unit  of  keq  is  dependent  on  n  itself: 

r  \  pressure  units] 

[k  units ]  =  — - j  (48) 

[sinkage  units]" 


To  mitigate  this  difficulty,  Reece’s  revision  of  Bekker’s  equation  (Equation  (16))  was  implemented  so  that  the  unit 
of  keq  is  always  in  pressure  units  and  does  not  depend  on  n  .  Applying  Equation  (16)  (instead  of  Equation  (15)) 

to  solve  for  the  parameters  resulted  in  very  little  n-keq  ’  correlation  for  our  data,  as  shown  in  Figure  1 1 .  For  each 
plate  size,  keq  ’  demonstrated  no  visible  trend  with  n  . 
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Figure  11:  Relationship  between  n  and  keq  \  Reece's  modification. 

For  each  experimental  run,  Equation  (16)  was  fit  to  a  plot  of  the  pressure  as  a  function  of  sinkage.  The  equivalent 
pres  sure- sinkage  constant,  keq  ,  and  exponent  of  sinkage  to  plate  width,  n  ,  found  in  Reece’s  equation  were 

determined  by  minimizing  an  error  function  based  on  the  least-squares  method  (analogous  to  method  3,  described  in 
Section  4. 1.2.1).  The  mean  and  standard  deviation  were  calculated  to  quantify  the  experimental  variability  of  each 
soil  parameter.  The  parameters  for  each  data  set  are  recorded  in  Table  3. 

Table  3:  Pressure-sinhage  parameters  and  the  mean  and  standard  deviation  for  each  plate  size,  Reece's  equation. 

Values  from  experimental  runs  are  italicized. 


Plate  Width, 
b  (cm) 

3 

5 

7 

Test# 

n 

keq/  [MPa] 

n 

keq/  [MPa] 

n 

keq,/  [MPa] 

1 

0.938 

0.070 

1.208 

0.140 

1.082 

0.231 

2 

0.876 

0.067 

0.784 

0.134 

1.114 

0.266 

3 

1.289 

0.059 

0.774 

0.149 

0.976 

0.252 

4 

1.017 

0.065 

0.839 

0.133 

0.996 

0.250 

5 

1.162 

0.066 

0.773 

0.136 

1.053 

0.210 

6 

0.928 

0.068 

0.983 

0.160 

0.960 

0.217 

7 

0.980 

0.066 

0.775 

0.132 

1.154 

0.236 

8 

0.958 

0.070 

0.848 

0.125 

1.028 

0.222 

9 

1.434 

0.048 

0.716 

0.139 

1.066 

0.270 

10 

0.930 

0.073 

1.053 

0.137 

1.224 

0.327 

11 

0.979 

0.061 

0.874 

0.157 

1.085 

0.288 

12 

1.002 

0.067 

0.938 

0.162 

1.102 

0.245 

13 

0.926 

0.070 

1.130 

0.119 

1.107 

0.254 

14 

1.079 

0.057 

0.797 

0.147 

1.252 

0.262 

15 

1.061 

0.057 

0.851 

0.851 

1.440 

1.440 

Mean 

1.037 

0.064 

0.889 

0.141 

1.109 

0.255 

Std.  Dev. 

0.152 

0.007 

0.145 

0.013 

0.123 

0.032 
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The  standard  deviation  of  the  Reece  parameters  were  propagated  to  uncertainty  in  pressure  by  applying  the 
statistical  techniques  described  previously  Section  4.2.3.  The  uncertainty  of  Reece’s  model  for  the  pressure,  cr ,  at  a 
given  sinkage,  z  ,  can  be  obtained  as: 


where 


(49) 


(50) 


and  the  uncertainties  of  the  independent  parameters  u  .  and  u  are  made  equivalent  to  the  standard  deviation  of  the 

keq 

pres  sure- sinkage  parameters  in  Table  3. 

Figure  12  -  Figure  14  show  the  uncertainty  estimation  of  the  pressure-sinkage  model  for  plate  widths  of  3,  5,  and  7 
cm,  respectively,  overlaid  with  the  mean  and  standard  deviation  of  the  experimental  pressure  data. 

As  can  be  seen  in  Figure  12  -  Figure  14,  although  Reece’s  model  does  not  precisely  capture  the  pressure-sinkage 
behavior,  the  model  uncertainty  provides  a  close  estimate  of  the  experimental  standard  deviation. 


Figure  12:  Uncertainty  estimation  for  Figure  13:  Uncertainty  estimation  for 
pressure-sinkage  test  for  3  cm  plate.  pressure-sinkage  test  for  5  cm  plate. 


Figure  14:  Uncertainty  estimation  for 
pressure-sinkage  test  for  7  cm  plate. 


In  contrast,  if  the  uncertainty  were  to  be  determined  by  the  standard  deviations  of  Bekker’s  parameters,  the  bounds 
would  be  much  wider  than  the  experimental  deviations,  as  illustrated  in  Figure  15  -  Figure  17. 
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parameters.  parameters.  parameters. 


The  mean  values  of  Reece’s  model  parameters  for  the  5  cm  plate  (keq>5  ’  =  0.141 ,  n  =  0.889)  were  used  to  estimate 
the  mean  pres  sure- sinkage  curve  for  the  wheel  performance  analyses  in  this  study.  This  plate  size  was  selected 
because  its  contact  patch  area  most  closely  resembles  that  of  the  single  wheel  used  in  this  study. 

5.1.2  Direct  Shear  Test 

A  total  of  12  direct  shear  tests  were  performed  to  analyze  the  shear  characteristics  of  the  dry  soil.  The  tests  were 
repeated  three  times  for  four  different  normal  stresses  each  (2.08,  2.86,  5.33,  and  17.83  kPa).  For  each  experimental 
run  that  was  carried  out  for  the  direct  shear  test,  Equation  (23)  was  fit  to  a  plot  of  the  shear  stress  as  a  function  of 
shear  displacement.  The  residual  shear  stress,  zres ,  and  shear  modulus,  K  ,  found  in  Equation  (23)  were  determined 

by  minimizing  an  error  function  based  on  the  least  squares  method.  The  parameters  for  each  data  set  are  recorded  in 
Table  4.  The  mean  and  standard  deviation  are  calculated  to  quantify  the  variability  of  each  soil  parameter. 

Table  4:  Soil  parameters  based  on  shear-displacement  curve  fits  along  with  the  mean  and  standard  deviation  for  each  parameter. 

Values  for  experimental  runs  are  italicized. 


Normal  Stress 
[kPa] 

2.08 

2.86 

5.33 

17.83 

Density  [g/cm3] 

1.5  (Loose) 

1.7  (Dense) 

1.7  (Dense) 

1.7  (Dense) 

*max  [Pal 

K  [m] 

*max  tPa] 

K  [m] 

*miuc  [Pa] 

K  [m] 

*max  [Pa] 

K  [m] 

Test  1 

2431.590 

5.010E-04 

2763.430 

1.021E-04 

4252.342 

2.073 E-04 

11364.000 

1.436E-04 

Test  2 

2468.297 

6.126E-04 

3189.900 

1.303 E-04 

3425.338 

2. 096 E-04 

13453.684 

1.297E-04 

Test  3 

2541.600 

6.300E-04 

1952.265 

2.260E-04 

3634.842 

2.528E-04 

12595.687 

1.1 00  E-04 

Mean 

2480.496 

5.812E-04 

2635.198 

1.528E-04 

3770.841 

2.232E-04 

12471.124 

1.277E-04 

Std.  Dev. 

56.010 

7.000E-05 

628.703 

6.494E-05 

429.948 

2.558E-05 

1050.396 

1.689E-05 

The  results  of  the  analysis  in  Table  4  were  used  to  determine  the  cohesion  and  angle  of  friction  of  the  soil.  The 
maximum  shear  stress  and  applied  normal  stress  are  related  using  Equation  (27),  and  are  shown  in  Figure  18. 
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Figure  18:  The  cohesion  and  angle  of  friction  of  the  soil  are  determined  using  the  Mohr-Coulomb  failure  criteria.  The  cohesion 

is  the  intercept  of  the  plot  and  the  angle  of  friction  is  the  slope. 

Based  on  this  linear  curve  fit  and  Equation  (27),  the  cohesion  of  the  material  is  equal  to  714.971  Pa  and  the  angle  of 
internal  friction  is  equal  to  33.083  degrees.  Using  the  techniques  described  in  Section  4.2.3  and  Equation  (28),  the 
uncertainty  is  415.491  Pa  and  2.875  degrees  for  cohesion  and  internal  friction  angle,  respectively. 

To  demonstrate  the  propagation  of  the  variability  of  the  experimental  parameters  to  the  uncertainty  of  Janosi  and 
Hanamoto’s  model,  the  techniques  in  Section  4.2.3  were  applied  using  the  data  in  Table  4.  Figure  19  -  Figure  22 
compare  the  uncertainty  estimation  of  the  direct  shear  tests  for  varying  normal  stresses  to  the  mean  and  standard 
deviation  of  the  experimental  data.  Although  Janosi  and  Hanamoto’s  model  does  not  precisely  capture  the  shear- 
displacement  behavior  for  dense  soil,  the  cohesion  and  angle  of  friction  have  been  shown  to  be  independent  of 
density  when  using  the  residual  shear  stress  [10].  It  is  important  to  note  that  the  shear  modulus  is  dependent  on 
density.  Janosi  and  Hanamoto’s  model  was  used  to  obtain  the  mean  value  of  the  shear-displacement  estimation.  The 
uncertainty  of  Equation  (23)  for  the  shear  stress,  r  ,  at  a  given  shear  displacement,  j  ,  is  approximated  as: 


2  2 

U  =  U 


dr 

Sr 


+  uv 


dr 

dK 


(51) 


where 


dr 


dr 


res 


dr 

dK 


(52) 


and  the  uncertainties,  uT  and  uK  ,  of  the  independent  parameters  are  determined  by  the  standard  deviation  of  the 
parameters  determined  by  the  direct  shear  tests  in  Table  4. 
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Figure  19:  Model  uncertainty  estimation  and  experimental 
standard  deviation  for  a  normal  stress  of  2.08  kPa  on  loose 
soil. 


Figure  21:  Model  uncertainty  estimation  and  experimental 
standard  deviation  for  a  normal  stress  of  5.33  kPa  on  dense 
soil. 


Figure  20:  Model  uncertainty  estimation  and  experimental 
standard  deviation  for  a  normal  stress  of  2.86  kPa  on  dense 
soil. 


Figure  22:  Model  uncertainty  estimation  and  experimental 
standard  deviation  for  a  normal  stress  of  17.83  kPa  on  dense 
soil. 


5.2  Wheel  Performance  Uncertainty 

Force  sensors  at  five  locations  across  the  width  of  the  wheel  (Figure  6)  measured  the  normal  and  tangential  stress 
distributions  (Figure  23  -  Figure  28).  Similarly,  the  drawbar  pull,  torque,  and  wheel  sinkage  at  five  different  slip 
ratios  ranging  from  -70%  to  70%,  and  two  different  normal  loads  of  70  N  and  135  N,  shown  in  Figure  29  and  Figure 
30,  were  measured  using  the  single-wheel  test  bed  in  Figure  4. 

The  uncertainty  estimation  method  described  in  Section  4.2.3  is  used  along  with  the  soil  parameters  in  Table  5  and 
wheel  properties  in  Table  6  to  compare  the  wheel  performance  model  to  the  experimental  values  of  wheel 
performance.  The  normal  and  tangential  stresses  at  the  wheel-soil  interface  were  calculated  from  Equations  (16)  and 
(30),  respectively,  and  plotted  in  Figure  23  -  Figure  28.  The  drawbar  pull,  torque,  and  wheel  sinkage  at  three 
different  slip  ratios  ranging  from  10%  to  70%,  and  two  different  normal  loads  of  70  N  and  135  N,  shown  in  Figure 
29  and  Figure  30,  were  calculated  using  Equations  (37),  (38),  and  (39). 

Table  5:  Soil  parameters. 


Parameter 

Value 

Uncertainty 

Units 

Pressure-sinkage  constant,  k  b=5cm 

1.41xl05 

1.30xl04 

Pa 
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Exponent  of  (sinkage\width),  nb=5cm 

0.889 

0.145 

- 

Cohesion,  c 

714.971 

415.491 

Pa 

Angle  of  internal  friction,  (p 

33.083 

2.875 

degrees 

Shear  modulus,  K 

5.571xl0‘4 

7.525xl0‘5 

m 

Coefficient  for  maximum  stress,  c1 

0.5 

0.2 

- 

Coefficient  for  maximum  stress,  c2 

0.5 

0.2 

- 

Table  6:  Wheel  properties. 


Parameter 

Value 

Uncertainty 

Units 

Normal  force,  Fz 

70,  135 

5% 

N 

Wheel  radius,  r 

0.13 

0 

m 

Wheel  width,  b 

0.16 

0 

m 

The  uncertainty  in  the  wheel  performance  can  be  calculated  using  the  techniques  described  in  Section  4.2.3.  The 
model  uncertainties  of  the  normal  and  tangential  stresses  for  W  =  100  N  are  plotted  against  the  experimental  data  in 
Figure  23  -  Figure  28.  The  uncertainty  bars  are  determined  numerically  [16]  using  a  finite  differencing  technique  to 
obtain  the  sensitivity  of  the  wheel  performance  with  respect  to  each  variable.  Note  that  the  model  uncertainty 
envelope  does  not  sufficiently  capture  the  experimental  data. 


Figure  23:  Normal  stress  vs.  contact 
angle  for  10%  slip. 


Figure  26:Tangential  stress  vs.  contact 
angle  for  10%  slip. 


Figure  24:  Normal  stress  vs.  contact 
angle  for  30%  slip. 


Figure  27:  Tangential  stress  vs.  contact 
angle  for  30%  slip. 


Figure  25:  Normal  stress  vs.  contact 
angle  for  70%  slip. 


Figure  28:  Tangential  stress  vs.  contact 
angle  for  70%  slip. 


Finally,  the  uncertainty  of  the  soil  parameters  is  propagated  to  the  drawbar  pull,  torque,  and  sinkage  of  the  wheel,  as 
shown  in  Figure  29  and  Figure  30.  The  uncertainty  bars  are  determined  numerically  [16]  using  a  finite  differencing 
technique  to  obtain  the  sensitivity  of  the  wheel  performance  with  respect  to  each  variable  and  are  plotted  along  with 
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the  measurements  based  off  the  individual  experimental  runs.  Note  that  the  model  uncertainty  envelope  does  not 
sufficiently  capture  the  experimental  data,  particularly  in  torque  and  sinkage. 


Figure  29:  The  uncertainty  of  the  drawbar  pull,  torque,  and  sinkage  of  the  wheel  for  a  plate  size  of  5  cm  x  16  cm  and  vertical 

force  of  70N. 
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Figure  30:  The  uncertainty  of  the  drawbar  pull,  torque,  and  sinkage  of  the  wheel  for  a  plate  size  of  5  cm  x  16  cm  and  vertical 

force  of  13 5N. 


6  Summary 


6. 1  Soil  Parameter  Variability 

The  first  objective  of  this  paper  was  to  characterize  the  variability  of  measured  soil  parameters  given  experimental 
data.  To  determine  the  soil  parameters  associated  with  normal  loads,  the  pressure  and  sinkage  were  measured  with  a 
plate  that  penetrated  the  soil  under  controlled  test  conditions.  To  determine  the  soil  parameters  associated  with 
shearing  loads,  a  standard  direct  shear  test  was  performed  under  various  normal  loads.  The  individual  soil 
parameters  were  obtained  from  the  pres  sure- sinkage  tests  by  fitting  the  dimensionless  form  of  Reece’s  equation  to 
each  experimental  data  set.  Similarly,  Janosi  and  Hanamoto’s  equation  for  shearing  was  fit  to  the  direct  shear  test  to 
determine  the  soil  parameters  for  each  experimental  data  set.  The  mean  and  standard  deviation  were  calculated  to 
quantify  the  variability  of  each  soil  parameter.  This  experimental  variability  was  used  to  calculate  the  uncertainty  of 
the  soil  stress.  Based  on  Figure  12  -  Figure  14  and  Figure  19  -  Figure  22,  the  model  uncertainties  show  good 
agreement  with  the  actual  standard  deviations  of  the  experimental  measurements.  Since  the  model  parameters  were 
estimated  from  the  same  experimental  data,  this  uncertainty  correlation  validates  the  uncertainty  propagation  theory 
employed  in  this  paper. 


6.2  Wheel  Performance  Uncertainty 

The  second  objective  of  this  contribution  was  to  quantify  the  uncertainty  in  wheel  performance.  Such  an 
investigation  requires  that  the  wheel  can  be  accurately  modeled.  The  sinkage  of  the  rigid  wheel  model,  based  on  the 
work  of  Wong  and  Reece,  was  determined  by  relating  the  weight  of  the  wheel  and  the  stresses  acting  on  it  in  the 
vertical  direction.  Given  the  uncertainties  of  the  individual  soil  parameters,  it  was  possible  to  determine  the  overall 
uncertainty  of  the  wheel  performance  (drawbar  pull,  torque,  and  sinkage)  using  the  uncertainty  propagation 
technique  described  in  Section  4.2.3.  The  uncertainties  of  the  normal  and  tangential  stresses  were  plotted  against  the 
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experimental  data  for  a  single  wheel  traversing  over  the  soil  in  Figure  23  -  Figure  28.  Additionally,  the  uncertainty 
of  the  soil  parameters  was  propagated  to  the  drawbar  pull,  torque,  and  sinkage  of  the  wheel,  as  shown  in  Figure  29 
and  Figure  30.  The  wheel  performance  model  fails  to  capture  much  of  the  experimental  stress  data  and  has  little 
overlap  with  the  drawbar  pull,  torque,  and  sinkage.  It  is  evident  from  the  results  that  the  classical  terramechanics 
model  used  in  this  paper  is  not  suitable  for  the  analysis  of  lightweight  vehicles  that  exhibit  relatively  low  ground 
pressures  (i.e.  15-35  kPa). 

Several  modifications  could  be  made  improve  the  accuracy  of  this  model  for  lightweight  vehicles.  For  example,  the 
flat  plate  assumption  of  pressure-sinkage  models  has  been  shown  to  be  inaccurate  for  lightweight  vehicles  and 
would  likely  benefit  from  revisions  according  to  the  recent  results  for  diameter-dependent  models  reported  in  [17]. 
Further,  the  shear  characteristics  of  the  soil,  shown  in  Figure  19  -  Figure  22,  are  not  adequately  captured  and  could 
be  more  accurately  estimated  from  making  a  modification  to  Equation  (23)  such  as  in  [18]  to  account  for  the 
nonlinear  “hump”  that  occurs  in  dense  soil.  Lastly,  based  on  the  size  of  the  vehicle  and  the  nature  of  the  terrain, 
alternate  methods  could  be  employed  to  capture  the  behavior  of  each  individual  soil  particle  [19]. 


7  Appendix  A:  Matlab  Code  for  WheelSoil  Class 


The  WheelSoil  class  encapsulates  the  Bekker-Wong  Theory  for  vehicle-soil  interaction.  A  WheelSoil  object 
contains  the  necessary  variables  and  functions  to  compute  wheel  stress  and  performance,  given  the  parameters  of  the 
soil  and  wheel. 
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classdef  wheel soil  <  handle 

%  Encapsulates  Bekker-Wong  Theory  for  Vehicle-Wheel  Interaction 
%  input:  Soil  parameters ,  Output:  Wheel  performance 


properties 

%  properties  with  temporary  values 


phi_degree  ■  19,735; 
phi  ■  0; 
c  =  293,123; 

K  =  0,0006350 9; 
k_l  *  0 , 141e6 ; 
k_2  ■  0; 
n  -  0,39; 
gamma  ■  0,0575; 


%  angle  of  internal  shearing  resistance 
%  cohesion.  Pa 

%  shear  deformation  modulus,  m 
%  pressure  sinkage  modulus  1 
%  pressure  sinkage  modulus  2 
%  exponent  of  sinkage  to  width  ratio 
%  density,  N/irr3 


%  Wheel  properties 
r  *  ,13; 
b  -  ,16; 

W  *  135; 


%  radius ,  m 
%  width,  m 

%  vertical  axle  load,  N 


%  Coefficients  for  the  relative  position  of  max,  radial  stress 
c_l  =  0,5;  %  0,43; 

c_2  a  0,5;  %  0,32; 

slip  *0,7;  %  slip 


%  important  angles  along  the  wheel 
theta_l  ■  0,1; 

theta_2  =0;  %  assume  no  t 

theta_m  a  Q; 


z_0  ■  0; 

end 


%  maximum  sinkage,  m 
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methods 

%  Constructor  for  the  Wheel Soil  class  (contains  the  variables) 
function  TS  ■  Wheel£oil(±,ki&q,n,c,pihi,K,cl,c2,Fz,r,b,jfcheta2) 

TS ,phi_dcgrcc  -  phi?  %  angle  of  internal  shearing  resistance 
TS,phi  *  0? 

TS,c  *  c?  %  cohesion.  Fa 

TS,K  =  K?  %  shear  deformation  modulus,  m 


T5,k_l  *  keq? 
TS,k_2  *  0? 
TS,n  =  n? 


%  pressure  sink age  modulus  1 
%  pressure  sinkage  modulus  2 
%  exponent  of  sinkage  to  width  ratio 


TS, gamma  »  0,043? 


%  density,  N/m~3  (doesn't  matter) 


%  Wheel  properties 
TS,r  ■  r? 

TS,b  -  b? 

TS,W  *  Fz? 


%  radius,  m 
%  width,  m 

%  vertical  axle  load,  N 


%  Coefficients  for  the  relative  position  of  max,  radial  stress 

TS ,C_1  ■  Cl? 

TS ,c_2  ■  c2? 

TS,slip  =  i?  %  slip 

TS,theta_l  *  0,1; 

TS,theta_2  =  theta2*pi/130 ?  %  assume  no  t 

TS , t  het  a_m  =  0? 

TS , z_0  ■  0? 

TS,slip  h  i+ 

TS,phi  ■  TS,phi_degree*pi/,13Q? 

end 


%  Get  the  normal  stress  along  the  front  of  the  wheel  for  a  given 
%  angle 

function  sigma_l  ■  getSigmal(TS, theta) 

z  ■  ( cos ( t  het  a ) - cos ( TS ,  t  het  a_l ) ) *T5 , r ? 
b_plate  =  0,05? 

s i gma_l  ■  TS , k_l , * ( z , /b_pl ate ) . ~TS , n  ?  %  k_l  is  k_eq  here 

end 


%  Get  the  normal  stress  along  the  back  of  the  wheel  for  a  given 
%  angle 

function  sigma_2  ■  get Sigma2(TS, theta) 

z  ■  (cos(TS,theta_l-( (theta-T5,theta_2)/ . , , 

(TS , t.heta_m-T5 , theta_2 ) ) * (TS , theta_l-TS , thetajn) ) .  r . 
-cos (TS , theta_l ) ) *TS , r ? 
b_plate  ■  0,05? 

s i gma_2  ■  TS , k_l , * ( z , /b_pl ate ) - ^TS , n  ?  %  k_l  is  k_eq  here 

end 
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%  Get  the  normal  stress  at  any  angle  on  the  wheel 
function  sigma  =  get5igma{T5, theta) 
sigma  *  zeros (length (theta) , 1) j 
for  i*li Length (theta) 

if (theta (i)>=TS.theta_in  &&  theta(i)<*TS. theta_l) 
sigma (i)  =  get5igmal(T5, theta(i) ) ; 
el s  eif ( thet  a { i ) > =T5 , t het a_2  &  i  t het a ( i ) <TS , thet a_m ) 
sigma(i)  =  get5igma2(TS,theta(i) ) j 

else 

sigma{i)  ■  0; 

end 

end 

sigma  ■  sigma' ; 

end 


%  Get  the  shear  stress  at  any  angle  on  the  wheel 
function  tau  *  getTau(TS, theta) 
tau  ■  zeros {length( theta) ,1) ; 
f  or  i®  1 : l  engt  h ( t  het  a ) 

if ( thet  a ( i ) >*TS  *  theta_2  &&  theta ( i ) <»TE . thet  a_l ) 
sigma  ■  TS.getSigma(theta(i) )  *f 

j  ■  ( (TS,t.heta_l-theta(i)  )-(l-TS.slip)*(sin(TS,theta_l) . , , 
“Sin( theta(i) ) ) )*T5,r; 

tau(i)  ■  (T5,c+sigma*tan{T5.phi) )*(l-exp(-j/TS,K) ) f 

else 

tau(i)  a  0; 

end 

end 

tau  =  tau1  ; 

end 

%  Error  function,  including  derivative  (used  for  optimization) 
function  | error , derrordthet ac ]  =  er rorFunc t ion ( TE  , thet a_c ) 
error  ■  TS , getE  r ror ( thet  a_c ) ; 
if  (nargout  >  1) 

%  use  complex  differentiation 

derrordt  het  ac  *  imag { TS , getE  rror ( t  het  a_c + sqrt ( - 1 ) * 1 e- 1 0 ) ) / le- 1 0 

end 

end 


%  Calculate  the  normalized  error  in  weight  (used  for  optimization) 
function  error  ■  getError ( TS , thet a_c ) 

W_prime  =  ts , get Wei ght ( thet  a_c ) ; 
error  »  (W_prime-TS,W) j 

end 
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%  Get  the  weight  for  a  given  entrance  angle 
function  wcheck  *  getweight (TS, theta_c) 
radius  ■  TS.r; 
width  ■  TS.b; 

TS , thet  a_l  =  thet  a_c  ; 

T5.theta_iri  ■  (T5.c_l+T5.c_2*T5.slip)*theta_c; 

Wc  hec  k  ■  i  ntegr  al  ( I  ( thet  a )  ( TS .  ge  t  S  i  gma  { t  het  a ) .  *cos  ( thet  a ) )  ,  „ 
TS .  thet  a_2 ,  TS .  thet  a_l ) ; 

Wc hec  k  ■  We hec  k  +  i  nt  egr  al ( I ( t  het  a ) ( TS . getTau ( t  het  a ) . * . . . 

sin( theta) ) ,T3. theta_2,T5. theta_l) ; 
wcheck  ■  wcheck.* radius.* width; 

end 

%  Calculate  the  wheel  performance  for  a  given  set  of  soil 
%  parameters 

function  [H,R,D,T,z]  ■  calculateWheelPerformance(TB , slip,keqfn,c 
phi, K, cl, c2,Fz,r,b, thet a2, guess) 

TS.slip  ®  slip; 

T5.k_l  =  keq; 

TS.n  =  n; 

TS .c  ■  c; 

TS.phi  ■  phi; 
ts .K  *  k; 

TS.c_l  ■  cl; 

T5.c_2  ■  c2; 

TS .W  =  Fz; 

TS.r  =  r; 

TS.b  ■  b; 

TS.theta_2  ®  theta2; 
radius  ■  r; 
width  »  b; 

%  Calculate  the  entrance  angle  (assuming  an  exit  angle) 
options  =  optimset ( 1  Jacobian 1 , 'on' , 1  Display 1 , 1  of f 1 ) ; 
thet  a_c  h  f  solve ( ITS . er rorFunc  t ion , gues  s , opt ions ) ; 

T5.theta_!  ■■  theta_c; 

TS.theta_m  -  (TS.c_l+TS.c_2*TS.slip)*TS.theta_l; 

TS . z_0  ■  ( 1-cos (TS. theta_l ) )*TS.r; 

%  Calculate  the  total  thrust 

H  a  int  egr  al ( $ ( t  het  a ) ( TS . getTau ( t  het  a ) . *  cos ( t  het  a ) ) ,  . . . 

TS . theta_2 ,TB . theta_l ) ; 

H  =  H.*radius.*width; 

%  Calculate  the  total  motion  resistance 

R  *  integral (|( theta) (TS.getSigma( theta) .*sin( theta) ) ,  r  r  , 

TS . theta_2,TS . theta_l ) ; 

R  B  R.*radius.*width; 
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%  Calculate  the  total  drawbar  pull 
D  ■  H-R? 

%  Calculate  the  total  torque 

T  ■  integral (@( theta) (T5,getTau( theta) ) ,TS,theta_2,TS,theta_l) ? 
T  H  T, *( radius” 2 * width ) ? 

^Calculate  the  maximum  sinkage 
z  ■  (1-cos (TS, theta_l) ) ,*radius? 

end 

end 

end 


8  Appendix  B:  Matlab  Code  for  Example  Program 


The  following  example  program  uses  the  Wheel  Soil  class  to  calculate  the  wheel  performance  and  uncertainty  for 
the  soil  and  wheel  described  in  this  paper.  Given  the  wheel  slip  and  wheel  weight,  the  program  will  return  the  values 
and  uncertainty  for  drawbar  pull,  torque,  and  sinkage. 


f unc t ion  [ D , T , z , D_unc , T_unc , z_unc ]  =  Wheel soilE xampl ©Driver (slip, wei ght ) 
%  Example  code  for  using  the  Wheel Soil  class 
%  Input;  wheel  slip 

%  weight  of  the  wheel 

%  Output:  drawbar  pullr  D 

%  torque ,  T 

%  maximum  sinkage t  z 

%  drawbar  pull  uncertainty ,  D_unc 

%  torque  uncertainty ,  T_unc 

%  sinkage  uncertainty ,  z_unc 

%  Wheel /soil  Parameters 


phi  *  19,735? 

% 

angle  of  friction,  degrees 

c  ■  293,123? 

% 

cohesion 

K  =  0,00063509? 

% 

shear  modulus 

k_eq  =  0.141e6? 

% 

pres sure- sinkage  modulus 

n  ■  0,09? 

% 

exponent  of  sinkage  to  width 

theta_2  =  -5? 

% 

exit  angle  (assumed) 

r  =  ,13? 

% 

wheel  radius 

b  ■  ,16? 

% 

wheel  width 

w  *  weight? 

% 

wheel  weight 

%  Coefficients  for  determining  the  relative  position  of  max,  radial  stress 
c_l  «  0.5? 
c_2  ■  0*5? 

%slip  ■  0,7? 
guess  *  0,1? 
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%  soil  parameter  variability 
phi_unc  *  0,02229; 
c_unc  =  239,7; 

K_unc  *  0,0001858; 
ft_eq_une  ■  13000; 
n_unc  *  0,14; 
theta_2_unc  ■  0,03; 

W_unc  =  2,43; 
c_l_unc  =  0,2; 
c_2_unc  *  0,2; 

%  Create  the  Wheel Soil  object 

test  =  WheelSoil{ slip, k_eq,n,c, phi, K,c_l,c_2,W,r,b,theta_2) ; 

%  Calculate  the  uncertainty! 

%  Need  the  partial  derivatives  w.r.t.  soil  parameters 

%  Calculate  the  partial  derivatives  w.r.t,  pressure  sinkage  modulus 
[  H ,  R ,  D ,  T ,  z  ]  *  test.c ale ul ateWheelPerf ormance ( s 1 ip , k_eq+  sqr  t ( - 1 )  * le- 1 0 , r 
n, c  ,phi*pi/180  ,K, c_l  ,c_2  , W,  r , b,  theta_2*pi/18Q  ,  guess ) ; 
dDdk_eq  *  imag(D)/le-10; 
dTdk_eq  ■  imag{T)/le-10; 
dzdk_eq  =  imagi;z)/le-10; 

%  Calculate  the  partial  derivatives  w.r.t.  exponent  of  sihkage  to  width 
[  H ,  R ,  D ,  T ,  z  ]  *  test.c  ale  ul at eWheelP er  f ormanc e { s 1 ip  , k_eq  f  r  r  r 

n+sqrt ( -1 ) *le-10 ,c ,phi*pi/130 ,K,c_l ,c_2  ,W, r , b, theta_2*pi/180 ,  guess ) 
dDdn  ■  imag{D)/le-10; 
dTdn  =  imag{T)/le-10; 
dzdn  *  imag(z)/le-10; 

%  Calculate  the  partial  derivatives  w.r.t.  cohesion 
[H,R,D,T,z]  =  tes  t . c  ale  ul ateWheelPerf ormance ( s 1 ip  , k_eq  , n  ,  „  .  * 

c+sqrt (-1 ) *le-10 ,phi*pi/180 ,K,c_l , c_2  , W, r , b, theta_2*pi/180  , guess ) ; 
dode  ■  imag(D)/le-10; 
dTdc  *  imag{T)/le-10; 
dzde  ■  imag{z)/le-10; 

I  Calculate  the  partial  derivatives  w.r.t.  angle  of  friction 
[Ei,R,D,T,z]  «  test.c ale ulateWheelPerformanee( slip, k_cq,n,c,  -  .  . 

phi*pi/130+sqrt ( -1 ) *le-10 , K,c_l , c_2 ,W, r ,b, theta_2*pi/130  , guess ) ; 
dDdphi  a  imag(D) /le-10 ; 
dTdphi  =  imag{T) /le-10 ; 
dzdphi  ■  imag(z) /le-10; 
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%  Calculate  the  partial  derivatives  w.r.t.  shear  modulus 
[ H,Rf  D^T, z J  *  tes t . c ale ul at  eWheelPe r f ormanc e ( s lip , k_eq ,n,cf  .  .  . 

phi*pi/18Q ,K+sqrt ( -1 } *le-10 , cl , c_2  , W, r r hf tharta_2*pi/10Q  , guess ) 
dDdK  =  imag{B)/le-10; 
dTdK  ■  imag{T) /le-10; 
dzdK  ■  imag{z) /le-10; 

%  Calculate  the  partial  derivatives  w.r.t.  constant  1 
[H,RfDlrTfzj  =  tes  t ,  c  ale  ul  at  eWheelPer  f  ormanc  e  ( s  1  ip  ,  k_eq  ,  n  ,  c  ,  *  *  . 

phi*pi/180 c_l+sqrt ( -1 ) *le-10 ,c_2  r , b, theta_2*pi/18Q ,  guess ) 

dDdc_l  *  imag{D) /le-10; 
dTdc_l  ■  imag(T) /le-10; 
dzde_l  =  imag(z) /le-10; 

%  Calculate  the  partial  derivatives  w.r.t.  constant  2 
[  H ,  R ,  D ,  T ,  z  ]  *  tes  t , c ale  ul at eWheelPer f ormanc e ( s 1 ip , X_eq ,n,Gf , „ , 

phi*pi/180  , K, c_l , c_2+sqrt ( -1 ) * le-10  , W, r , b, theta_2*pi/180 ,  guess ) 
dDde_2  ■  bmag(D) /le-10; 
dTdc_2  =  imag{T) /le-10 ; 
dzde_2  «  imag(z) /le-10; 

%  Calculate  the  partial  derivatives  w.r.t.  wheel  weight 
[H,RiD,T,z]  »  tes  t ..  c  ale  ul  at  eWheelPerf  ormanc  e  ( s  1  ip  ,  k_eq  ,  n  ,  c  ,  .  .  . 

phi *pi/ 1 8  0  , K  , c_l , c_2 , w+  sqrt ( - 1 ) * le- 1 0  , r , b  r t het a_2  *pi/ 180,  guc  s  s ) 
dDdW  a  imag{D) /le-10; 
dTdW  a  imag{T) /le-10; 
dzdw  ■  imag{z) /le-10; 

%  Calculate  the  partial  derivatives  w.r.t.  exit  angle 
[  H ,  R ,  D ,  T ,  z  ]  =  t  es  t , e  ale  ul at eWheelPer f ormanc e ( s 1 ip  , k_eq  , n  , c  ,  *  *  r 

phi*pi/180  ,K, c_l ,c_2 ,W, r , b, theta_2*pi/180+sqrt ( -1 ) *le-10 ,  guess ) 
dDdtheta_2  *  imag(D) /le-10; 
dTdtheta_2  -  imag{T) /le-10; 
dzdtheta_2  ■  imag(z) /le-10; 

%  Calculate  the  total  uncertainty  in  wheel  performance 

B_unc  *  sqrt { dcdk_eq “ 2  * k_eq_unc " 2  +dDdn  ~ 2  *  n_unc “ 2  +dDdc 2  *c_unc  * 2  +  « . . 

dDdphi " 2  *phi_unc " 2  +dDdK " 2  *K_unc " 2  +dBde_l * 2  *c_l_unc * 2  + _ 

dDdc_2 ” 2  *e_2_unc 2  +dBdW " 2  *W_unc * 2  +dDdt  het  a_2  ”2*1 het  a_2_unc " 2 ) ; 
T_unc  =  sqrt { dTdk_eq '2* k_eq_unc ” 2  +dTdn " 2  * n_unc ” 2  +dTdc " 2  *c_une ” 2  + . , _ 
dTdphi " 2  *phi_unc  ~ 2  +dTdK  ~ 2  *K_unc  * 2  +dTdc_l " 2  *c_l_unc " 2  + . . . 
dTde_2  ”  2  *e_2_unc "  2  +dTdW  ”  2  *’h'_unc "  2  +dTdt  het  a_2  *  2  *  t  het  a_2_unc  ’  2 ) ; 
z_une  ■  sqrt { dzdk_eq “ 2 * k_eq_unc ” 2  +dzdn “ 2 * n_unc ” 2  +dzdc " 2 *e_unc "2+ . . . 
dzdphi  ”  2  *phi_unc  "  2  +dzdK ~  2  *K_unc  "  2  +dzdc_l "  2  *c_l_unc 2  + . , , 
dzdc_2 ” 2  *e_2_unc  * 2  +dzdw ” 2  *W_unc * 2  +dzdthet a_2  ~ 2  *  t het a_2_unc " 2 ) ; 

%  Calculate  the  wheel  performance 

[ H f R f D , T f z J  «  tes t - c ale ul at eWheelPerf ormanc e ( s 1 ip , k_eq fn,c, . . . 
phi*pi/180 c_l ,c_2 r theta_2*pi/180 , guess  j ; 


end 
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